% 本程序演示了一个简单的Verlet算法。
% Verlet算法是一种经典的数值求解牛顿第二定律的方法，原则上可以模拟任意N粒子相互作用的问题。
% F=ma => a=F/m => r'' = F/m
% 随后代入二阶导数的离散形式，即可得到基本形式的Verlet算法。
% 默认演示了一个限制性三体问题。你可以改变粒子的初始状态，以及粒子间作用力的表达式。
% 然而，本实现比较粗糙，运算效率、精度较低，在运行一定时间后由于误差积累结果将失真。
% CC-BY 4

clc
clear

global n
global dt

n = 3; %粒子个数
dt = 0.01;

%---Three body
a = 2;
v = sqrt(10/a);

points(1,2).x = 0; %x, y: 粒子初始位置; vx, vy：粒子初始速度; m: 粒子质量
points(1,2).y = 0;
points(1,2).m = 1;
points(1,2).vx = v*cosd(60);
points(1,2).vy = -v*sind(60);

points(2,2).x = a;
points(2,2).y = 0;
points(2,2).m = 1;
points(2,2).vx = v*cosd(60);
points(2,2).vy = v*sind(60);

points(3,2).x = a/2;
points(3,2).y = sqrt(3)/2*a;
points(3,2).m = 1;
points(3,2).vx = -v;
points(3,2).vy = 0;

for i = 1:n
  points(i,1).x = points(i,2).x - points(i,2).vx * dt;
  points(i,1).y = points(i,2).y - points(i,2).vy * dt;
  points(i,1).vx = points(i,2).vx;
  points(i,1).vy = points(i,2).vy;
  points(i,1).m = points(i,2).m;
end

function f = compute_f(point) %计算粒子间的作用力
  global n
  f_direction = zeros(2);

  for i = 1:n
    f(i,i).x = 0;
    f(i,i).y = 0;

    for j = 1:i-1

      f_direction(1) = point(j).x - point(i).x;
      f_direction(2) = point(j).y - point(i).y;

      dist = norm(f_direction);
      f_direction /= dist;

      f(i,j).x = 10*f_direction(1)/dist^2; %计算一个平行于二粒子连线的吸引力，以r^2衰减。电磁力与引力满足该形式。
      f(i,j).y = 10*f_direction(2)/dist^2;

      f(j,i).x = - f(i,j).x; %牛顿第三定律：粒子i对j的力=粒子j对i的力
      f(j,i).y = - f(i,j).y;
    endfor
  endfor
end

t(1) = 0;
t(2) = dt;

imgind = 0;
for tick = 3:10000
  t(tick) = (tick-1)*dt;

  f = compute_f(points(:,tick-1));
  for i = 1:n
    fi.x = 0;
    fi.y = 0;
    for j = 1:n
      if i == j
        continue
      endif
      fi.x += f(i,j).x;
      fi.y += f(i,j).y;
    end

    points(i,tick).m = points(i,tick-1).m;
    points(i,tick).x = 2*points(i,tick-1).x - points(i,tick-2).x + fi.x/points(i).m*(dt)^2; %Verlet 算法
    points(i,tick).y = 2*points(i,tick-1).y - points(i,tick-2).y + fi.y/points(i).m*(dt)^2;

    points(i,tick).vx = (points(i,tick).x - points(i,tick-1).x) / (dt);
    points(i,tick).vy = (points(i,tick).y - points(i,tick-1).y) / (dt);
  end

  if mod(imgind,10) == 1
    clf;
    hold on
    axis equal
    axis([-10 10 -10 10])
    for i =1:n
     scatter(points(i,tick).x,points(i,tick).y);
    endfor
    drawnow
    pause(0.01);
  end
  imgind++;
end


